package scu.maqiang.eigenvalue;

import java.util.Arrays;

import scu.maqiang.fes.BCType;
import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES2Q41;
import scu.maqiang.mesh.Mesh2Q4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.EigenSSolver;
import scu.maqiang.numeric.MVO;
import scu.maqiang.numeric.SRMatrix;

public class LapEigenValue2DQ41 {

	public static void main(String[] args) {
		Mesh2Q4 mesh = new Mesh2Q4().square2D(80, 80);
		mesh.scale(Math.PI, Math.PI);
		
		FES2Q41 fs = new FES2Q41(mesh);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		SRMatrix M = new SRMatrix(fs.GetNdof());
		
		int N = 10;
		double[] D = new double[N];
		double[][] V = new double[N][fs.GetNdof()];
		double[] x = new double[fs.GetNdof()];
		Arrays.fill(x, 1.0);
//		fs.assembleHeatStiff(new double[] {1.0}, BVPType.COMMON, A);
////		A.sort();
////		System.out.println(A);
//		fs.assembleHeatMass(new double[] {1.0}, BVPType.CONSIST_MATRIX, M);
		
		ScalarFunc cf = (xyz, label, param) -> label == 10? 1.0: 1.0;
		fs.assembleStiff(new ScalarFunc[] {cf}, null, BVPType.COMMON, A);
//		A.sort();
//		System.out.println(A);
		ScalarFunc rhof = (xyz, label, param) -> label == 10? 1: 1;
		fs.assembleMass(new ScalarFunc[] {rhof}, null, BVPType.CONSIST_MASS, M);
		//M.sort();
		//System.out.println(M);
		//fs.applyBC_MBN(A, x, Direct.All, 0.0, 1, 2, 3, 4);
		fs.applyBCZero(A, x, Direct.X, BCType.RCE, 1, 2, 3, 4);
		//System.out.println(A);
		EigenSSolver solver = new EigenSSolver();
		//solver.inversePowerTran_PLU(A, M, fs.getnBand(), x, 2.0, N, D, V);
		solver.inversePowerTran_PCGSSOR(A, M, x, 2.0, N, D, V);
		System.out.println(MVO.toString(D));
		String outputFilePreFix = "EigenValue";
		for(int i = 0; i < N; i++) {
			mesh.toTecplot(outputFilePreFix + (i+1) + ".dat", V[i]);
		}
		mesh.toTecplot("Eigenvalues.dat", V);

	}

}
